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ABSTRACT 


A system  of  time  dependent  integral  equations  are  derived 
and  then  are  analytically  demonstrated  to  be  capable  of 
treating  scattering  by  a dielectric  interface.  A finite 
difference  method  is  demonstrated  to  be  capable  of  determining 
the  fields  scattered  by  an  obstacle  having  an  edge,  by  com- 
paring a numerical  solution  to  a canonical  solution  for 
scattering  by  a perfectly  conducting  wedge.  Both  methods  are 
applied  to  the  dielectric  platform  model  which  has  both  a 
dielectric  interface  as  well  as  an  edge.  The  results  obtained 
are  in  close  agreement  and  we  choose  to  generate  production 
data  for  ATLAS  I related  parameters  by  employing  the  finite 
difference  method;  however,  this  should  not  be  taken  as  an 
endorsement  that  this  method  is  always  preferable.  As  part 
of  the  investigation  we  identify  problems  for  which  either 
method  would  be  preferred.  Finally,  we  present  time  dependent 
plots  of  the  electric  field  at  points  in  the  working  volume 
that  show  the  amount  of  distortion  caused  by  the  platform. 
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SECTION  I 


INTRODUCTION  AND  SUMMARY 

f ' 

The  primary  issue  that  pervaded  this  investigation  was 
the  question  of  how  to  achieve  confidence  in  the  numerical 
data  that  would  be  generated.  The  approach  used  was  to 
compare  the  results  obtained  by  two  dissimilar  calculational 
procedures  after  each  was  demonstrated  to  be  capedsle  of 
yielding  results  known  to  be  valid  for  separate  canonical 
problems.  The  two  approaches  are  a coupled  system  of  time 
dependent  integral  equations  and  an  appropriate  finite 
difference  method. 

The  concern  about  the  validity  of  the  results  is  due  to  the 
fact  that  the  plate  model  of  the  ATLAS  I trestle  platform  requires 
the  proper  numerical  treatment  of  a dielectric  Interface  as 
well  as  an  edge.  The  vast  majority  of  numerical  time 
dependent  scattering  calculations  deal  with  scattering  by  a 
perfectly  conducting  obstacle  and  few  of  those  studies  focus 
on  the  effect  of  singularities  caused  in  the  solution  due  to 
an  edge.  We  are  not  aware  of  £uiy  previous  solution  in  the 
literature  of  the  coupled  system  of  Integral  equations  that 
we  derive,  euid  numerically  solve.  The  finite  difference 
method  that  we  employ  has  been  previously  used  by  Page  emd 
Peterson  (ref.  1)  for  a dielectric  interface;  however,  the 
procedure  they  employ  at  the  Interface  is  different  from  our 
procedure.  This  difference  in  the  procedures  does  not  appear 
to  have  serious  consequence  since  numerical  testing  Indicated 
that  both  procedures  yield  similar  results  for  sufficiently 
small  grid  step  size. 


Page,  W.  E.  and  D,  H.  Peterson,  A Numerical  Method  for 
Computing  the  Propagation  of  an  Electromagnetic  Pulse 
Guided  Over  a Material  Interface,  Sensor  and  Simulation 
Note  96,  Air  Force  Weapons  Ledsoratory,  1970. 
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The  canonical  problems  used  to  test  the  system  of 
Integral  equations  are  the  problem  of  scattering  by  a 
dielectric  half  space  and  the  problem  of  scattering  by  an 
infinite  dielectric  slab.  For  these  problems,  we  analytically 
solve  our  system  of  integral  equations  and  obtain  the  known 
solutions . 

Our  testing  of  the  finite  difference  method  was  considerably 
more  involved  emd  a consequence  of  this  testing  has  the 
potential  for  yielding  significant  side  benefits.  We  utilize 
the  known  canonical  solution  for  scattering  of  a plane  wave 
step  fvinction  by  a perfectly  conducting  wedge.  This  solution 
was  convolved  with  the  function  of  time  that  we  Intended  to 
choose  for  an  Incident  plane  wave  pulse  that  would  be  conven- 
ient for  us  to  treat  by  the  finite  difference  method.  The 
result  of  this  convolution  describes  the  scattering  of  a plane 
wave,  having  the  desiredsle  time  dependence,  by  the  perfectly 
conducting  wedge.  It  is  a simple  matter  to  test  that  the 
results  obtained  by  the  convolution  procedure  are  accurate 
to  any  prespecified  number  of  significant  figures.  The  test 
consists  of  increasing  the  number  of  points  in  the  convolution 
integration  procedure.  We  compared  the  results  obtained  by 
the  finite  difference  method  with  the  results  of  known  accuracy 
obtained  by  the  convolution  approach  eind  determined  that  the 
finite  difference  approach  could  also  yield  solutions  to  any 
prespecified  number  of  significant  figures  by  decreasing  the 
finite  difference  grid  step  size. 

The  side  benefit  of  this  testing  of  the  finite  difference 
method  was  that  we  were  able  to  show  that  the  convolution 
solution  emd  the  finite  difference  solution  were  still  in 
agreement  when  we  let  the  wedge  angle  approach  zero.  This 
demonstrated  that  a particular  finite  difference  method  was 
cap2Uale  of  determining  the  field  scattered  by  a particular 
perfectly  conducting  open  surface  (the  semi-infinite  half 
pleme) . This  identifies  an  area  of  investigation  that  has 
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the  potential  to  satisfy  a long  standing  need  in  the  area  of 

EMP  interaction  and  coupling.  The  question  of  whether  an  I 

appropriate  finite  difference  method  is  capable  of  determining 

the  fields  scattered  by  a nonplanar  open  surface  merits  a r ! 

thorough  Investigation.  This  capedslllty  is  necessary  to  j 

quantify  errors  introduced  by  the  many  approximations  | ■ 

currently  employed  to  calculate  the  currents  and  voltages  at  i i 

I ) 

the  inputs  of  subsystems  contained  within  metallic  | 

enclosures  (missiles,  aircraft,  tanks,  ships,  etc.). 

I i 

Returning  to  the  question  of  confidence  in  the  numerical  !;  | 

data  we  present  for  the  model  of  the  ATLAS  I trestle  platform,  we  I | 

have  explained  how  we  concluded  that  the  integral  equation  approach  | ] 

was  ceipable  of  treating  the  dielectric  interface  problem  and  » I 

our  finite  difference  method  was  capad>le  of  treating  the  I ? 


edge  problem.  Our  final  test  was  to  apply  both  methods  to  | i 

the  dielectric  plate,  of  finite  extent  and  thickness,  which 

has  both  an  edge  and  a dielectric  Interface.  The  results  \ 

obtained  by  the  two  methods  were  in  agreement  auid  we  chose  to 

use  the  finite  difference  method  to  generate  the  production  | 

runs  that  had  parauneters  chosen  to  study  the  effect  of  i 

trestle's  platform.  For  this  particular  problem,  the  finite  | 

difference  method  was  chosen  due  to  computer  memory  consider-  j 

at ions;  however,  this  should  not  be  taken  as  an  endorsement 

that  this  method  is  always  preferable  to  use.  As  part  | 

our  investigation  we  have  identified  problems  for  which  | 

either  method  would  be  preferred. 

The  results  of  our  production  runs  show  that,  according 
to  our  simplified  model  of  the  platform,  the  fields  in 
TRESTLE'S  working  volume  are  clearly  distorted  by  the  platform. 

This  distortion  occurs  to  the  pulse  shape  as  well  as  to  its 
aunplltude.  As  the  observation  point  is  chosen  further  in 
from  the  leading  edge  of  the  platform,  the  distortion  persists 
for  larger  distances  above  the  platform.  Our  present  study 
is  limited  to  a distance  from  the  leading  edge  that  corresponds 
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to  25  platform  thicknesses.  The  observation  dlst2mce  is 
fxindamentally  limited  by  our  two-dimensional  modeling  of 
a three-dimensional  platform.  The  deeper  in  and  higher  up 
we  choose  to  observe,  the  sooner  we  sense  the  effects  of 
the  sides  of  the  platform  that  are  not  Included  in  the  two- 
dimensional  model.  Even  with  this  limitation,  our  data  is 
applicable  for  times  longer  than  the  time  the  incident  field 
requires  to  achieve  its  maximum  cunplltude. 

We  view  the  amount  of  distortion  exhibited  by  our  model 
and  calculations  as  demonstrating  a need  for  further  investi- 
gations, both  theoretical  and  experimental,  in  order  to  assess 
and  assist  the  threat  relatability  of  tests  that  will  be  per- 
formed in  ATLAS  I.  These  investigations  should  include  a more 
detailed  model  of  the  entire  support  structure  as  well  as 
interactions  with  test  objects  and  other  portions  of  the 
simulator . 
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SECTION  II 

FORMULATION  OF  THE  INTEGRAL  EQUATION  METHOD 
• r 

In  this  section  we  derive  the  system  of  integral  equations 
for  the  problem  of  scattering  from  a dielectric  cylinder  of 
infinite  length  and  arbitrary  cross  section.  The  conflgiira- 
tlon  of  Interest  is  depicted  in  figures  la  and  lb.  The 
Incident  electromagnetic  field  is  given  by 


® » Eq  f(t  - 

_inc  _ „inc  ^ . „inc  ^ 

- • ®y  ”z  ®z 


E^  COS  3 

o 


E_  sin  B 


- zVc^), 


f(t  - z'/c) 


(z  - Zq)  cos  3 + (y  “ y^)  sin  3 


where  e„,  e„/  e_  are  unit  vectors,  3 is  the  angle  defined  in 
figures  la  and  lb,  is  the  free-space  characteristic 
impedance,  c^  is  the  speed  of  light  in  free  space  and  x^,  y^, 
are  the  coordinates  of  a point  on  the  surface  of  the 
cylinders,  that  is  swept  by  the  Incident  wavefront  at  t » 0 


The  dielectric  cylinder  is  homogeneous  with  a dielectric 
permittivity  etnd  a magnetic  permeability  equal  to  the 
vacuum  permeeUalllty  y^.  Our  problem  is  clearly  two-dimensional, 
i.e.,  all  the  physical  quantities  of  interest  are  Independent 
of  X.  Under  these  circumstances  one  can  show  that  the 
scattered  electric  field  will  only  have  an  x-component  %rtiereas 
the  scattered  magnetic  field  will  lie  entirely  in  the  yz 
plane.  Maxwell's  equations  for  the  total  field,  incident 
plus  scattered,  are  then  reduced  to 
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(2) 


t 

e 


TT 


'W 


9H 

■5F  " 


(3) 


where  e is  equal  to  outside  the  body  and  equal  to  Cj^  Inside 

the  body  and  J^(y,z,t)  is  the  source  of  the  incident  field 
located  far  away  from  the  scattering  volume.  (When  e ■> 
the  source  term  in  Equation  3 should  be  set  equal  to  zero.) 

The  boundary  conditions  across  the  surface  of  the  cylinder 
are:  and  H * s continuous,  i.e.,  total  tangential  electric 

and  magnetic  fields  should  be  continuous  (fig.  2) . If  the 
Incident  wavefront  has  a sharp  front,  i.e.,  the  fields  are 
nonzero  there,  then  at  t » 0 there  is  a discontinuity  of  the 
fields  across  the  boundary.  We  will  assume  that  f(u)  is  a 
smooth  function  of  u and  define  it  more  precisely  as  we  treat 
our  equations  numerically  in  subsequent  sections.  For  the 
derivation  of  the  system  of  integral  equations  we  need,  as  we 
shall  see,  continuity  of  E and  3E  /dn  where  n is  the  outward 
normal  on  the  surface  of  the  cylinder.  Continuity  of  9E^/dn 
is  Inferred  by  noting  that  (fig.  2) 


3E  3E  3E 

■Sir  * ^yz^x  : " " “Sy"  "y  "z 


3E  3E^ 

(“sin  0)  + “g"  (cos  0) 


3E 


3E. 


" ■§?  ^’'z^  " “5#  ^®y^ 


“^o  It  <«z  h * ^ *y^  “ "»'o  Tt  5 


(4) 
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where  the  penultimate  step  employed  equations  2.  Since  H • s 
is  continuous  across  the  boundary  for  all  times  we  understand 
that  so  are  d/9t(H  • s)  and  dE  /3n.  (We  exclude  pathological 

f * 

functions . ) 

Next  we  proceed  to  derive  the  system  of  integral  equations 
that  solves  our  two-dimensional  scattering  problem,  i.e.,  it 
' allows  the  calculation  of  the  scattered  fields  inside  and 
outside  the  dielectric  body.  We  begin  with  the  wave  equation 
satisfied  by  that  can  be  derived  by  manipulating  equations 
2 and  3:  ~ 

■ -T  - F(£',t')  (5) 

y ^ c^  at'V 

where  p'  * y'e  + z'e_, 

£.  •*  y z 

'?  = Ex  (£',t')  » total  electric  field 
-1/2 

c * Cq  ■ ^^0^0^  outside  the  cylinder  and 
c ■ Cj^  * inside  the  cylinder, 

P(£',t')  ■ p^Caj^/at)  outside  the  cylinder  and 
P{£',t')  * 0 inside  the  cylinder. 

Next,  we  Introduce  the  two-dimensional  Green's  function  G 
satisfying 

- «(£'  - £>  - « 

' <6) 

G(£',£;  t',t)  is  defined  in  an  infin^e  free-space  or  an 
infinite  dielectric  medium  depending  on  whether  c ■■  c^  or 
c c^.  The  solution  to  equation  6 is 
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G(£',£;  t',t) 


u[±C(t'  - t)  - |£'  - £|] 

[c^Cf  - t)^  - (£'  - 

(7) 


where  0(x)  *1  for  x > 0 and  U(x)  ■ 0 for  x < 0.  The  plus 
sign  in  equation  7 corresponds  to  the  retarded  solution,  i.e., 
an  observer  at  p*  senses  at  t'  a disturbance  caused  by  a 
source  at  £ fired  at  the  retarded  time  t * t'  - |£'  - £|/c. 

The  minus  sign  in  equation  7 corresponds  to  the  advanced  solu- 
tion, i.e.,  an  observer  at  £'  senses  at  t'  a disturbance  caused 
by  a source  at  £ fired  at  the  advanced  time  t * t'  + |£'  - £|/c. 
This  solution  violates  causality.  Notice,  however,  that  if 
we  switch  the  observation  and  source  space-time  points  the 
advanced  solution  of  equation  7 becomes  the  retarded  solution 
for  the  problem  of  a disturbance  observed  at  (£,t)  and  caused 
by  a source  at  £'  fired  at  t'  « t - (£'  - £l/c.  This  observation 
will  be  utilized  later  on  when  we  derive  our  integral  relationshps . 


OB 


'«'(£', f) 


« 


-¥(£', f)  ^G^(£'»P;  t’,t)]dA'dt'  - 'P^"°(£,t)  - '!'(£,t) 


(8) 


\i[^ere  C is  the  contour  shown  in  figure  2,  n*  is  the  outward 
normal  (fig.  2) , S„  is  the  region  bounded  by  the  circle  at 
infinity  and  C and  is  the  incident  electric  field. 

The  derivation  of  equation  8 assumes  that  the  contour 
integral  at  infinity  (resulting  from  Green's  identity)  has 
been  set  equal  to  zero.  The  reason  is  that  Y euid  dY/3n' 
in  the  integrals  can  be  replaced  by  the  scattered  fields 
(one  can  see  this  by  applying  equation  8 for  '?  ■ 
for  any  finite  t'  or  they  are  zero.  Thus  the  integration 
in  t*  is  over  one  instant  only  (t*  » +«)  and  it  can  be  shown 
to  have  zero  contribution.  The  t' -integration  in  the  second 
integral  in  equation  8 cam  be  performed  explicitly.  At  t'  * - 
the  scattered  fields  are  zero  throughout  region  amd  at 
t'  » 00  have  gone  to  zero  smoothly  to  assure  that  the  integral 
over  S is  zero.  Thus  equation  8 can  be  rewritten  as 


00 

'r(£,t)  - 'i'^®(£,t)  - [¥(£', f) 


3^  (£',£;  t',t) 


- t',t)  ^'P(£',t')]  ds'df.  (9) 

As  we  can  see  from  equation  9 the  scattered  field  at  £ and 
t is  due  to  contributions  from  points  £'  on  the  contour 
firing  at  t'.  Thus  t*  must  be  less  tham  t and  if  we  recall 
equation  7 we  understand  that  G (£*,£;  t',t)  must  be  tadcen 
with  the  minus  sign  in  front  of  c(t'  - t) , l.e.,  it  is  the 
advanced  solution  of  equation  6 with  c « c^. 

For  the  region  Inside  the  contour  C we  can  apply  a similar 
procedure  and  aurrive  at  the  following  equation 
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/fr 

-00  C 


■5nT®i^£'»£'  t',t) 


- (£’,£;  t',t)  ^'i'(£',t')]  ds'df 


(10) 


where  C is  the  seune  contoiir  as  in  equation  9,  n*  is  the  out- 
ward normal  and  the  Green's  function  is  the  advanced  solution 
of  equation  6 with  c « 

In  order  to  obtain  our  system  of  integral  equations 
we  let  the  observation  point  £ approach  the  contour  C.  If 
the  contour  C is  smooth  then  the  singularity  due  to  3G/3n' 
at  £'  » p results  in  a term  ± (1/2) ’f  {£, t)  (plus  for  the  inside 
and  minus  for  the  outside)  and  equations  9 and  10  give 


5 fo<£’‘> 


00 

,t)  - f 


3G, 


o 3n' 


- G 


-o.  c 


o 3n^  / 


ds'df 


j ^ d»'df  (u) 

where  we  have  used  the  subscripts  "o"  and  "1"  to  denote  the 
outside  and  inside  of  the  cylinder  respectively.  One  cam 
show  that  a principal  value  Integration  over  C (resulting 
from  the  limiting  process  £-'’£'  froa  the  outside)  is  not 
necessary  because  the  kernel  in  the  integral  is  not  singular 
as  £'  -►  £ when  both  £ and  £'  lie  on  the  contour.  (A  contour 
with  sharp  corners  is  discussed  at  the  end  of  this  section.) 

If  we  recall  the  continuity  of  ¥ and  3¥/3n  as  we  cross 
the  boundary  contour  we  understand  that  'J'q  “ = ¥(£,t)  and 

i 

1 

f 


3'F^/3n  » 3'Pj^/3n  = 3'?/3n  where  £ is  on  the  contour  C,  Thus 
equations  11  represent  a system  of  integral  equations  which 
allow  the  determination  of  f and  3T/3n  on  the  boundary  contour 
C.  Once  these  quantities  are  known  one  can  employ  equations 
9 and  10  to  determine  ¥ outside  the  body  euid  inside  the  body 
respectively.  Then  Maxwell's  equations  2 allow  the  determina- 
tion of  H everywhere. 

In  order  to  cast  our  system  of  equations  into  a form 
eunened>le  to  numerical  treatment  we  employ  the  explicit  form  of 
G given  by  equation  7 and  manipulate  the  resulting  integrals 
to  eliminate  the  apparent  singular  behavior  that  results  from 
the  differentiation  of  G.  We  have. 


3G 


7'G 


A 


n' 


V 


.R 


3G  £ 

^ R 


n* 


vdiere  R » |£'  - £| 


(12) 


3G  ^ C ) 6[c(t  - t*)  - Rl 

"Sr  ■ 5?  )"  o 1/2 


[c^(t  - t')2  - R^] 


[c2(t.  t')2-R2] 


U(c(t  - t')  - Rl 


(13) 


If  we  combine  equations  12  and  13  we  can  write 


y **2?  n VilTI 

LLc^(t  - t* 


)=  - r']' 


•^t'-t-R/C 


^ ico  [c^(t  - t')^  - R^] 


57y  a[c(t  - t')  - Rl  dt' 

(14) 
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Recalling  that 


/■ 


' dt' 


Ic^(t  - t')2  - r2] 


372 


t - t' 


2 5 29  1/2 

ric‘^(t  - t')-^  - 


we  can  rewrite  the  integral  on  the  right-hand  side  of 
equation  14  as 


00 

- Tif 


'P(£\t')  (£'  - £)  . n' 

-37J  u[c(t  - t’)  - R]  dt' 


[C^(t  - t')^  - r2] 


00 

“ ■ ^r/  (£'  - £)  - S'  u[c(t  - f)  - HI  a|r 


(; 


t - t' 


R^[c^(t  - t')^  - R^l 


T77 


dt' 


^ (t  - f)  'F(£',t')  utc(t  - f)  - Rl  (£' 
5? 


- £) 
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r[c^(t  - t*)^  - R^l 


00 


(t  - f)  (£'  - £)  • n' 


•^[c^(t  - t')2  - r2] 


|u [c (t  - t' ) 


Rl 

R] 


jat. 


- ct6Ic(t  - t')  - R] 

(f  - t)  (£'  - £)  • h* 


-CB  R 


22  22 
-^[c-^Ct  - f)^  - R^] 


Olc(t  - t')  - 


Rl  ^ df 


1 

I? 


▼(£'  - t') (£'  - £)  . n' 


Rlc^(t  - t')^  - R^l 


T7T 


t'  - t - R/c 


(15) 
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If  we  compare  equation  14  to  equation  15  we  see  that 


I 


L 


(f  - t)  (£'  - P)  • 

R^[c^(t  - t')^  - R^l 


“ 3(jr  I 

T7?  Ip- 


dt' 


(16) 


If  we  now  return  to  our  system  of  ecfuations  11,  we  can  use 
equation  16  to  rewrite  them  as 


§ 'F(£,t)  - 'i'^"^(£,t)  + ^ 

—00  C L 


) (f  - t) 


34' 

3t^ 


34* 

3rP 


1 U[c_(t  - f)  - R] 

■J  71 i Ipi 

-I  [c‘(t  - f)^  - r2 


dt' 


1 T(£,t)  - - ^ 


) (f  - t) 


U[c^(t  - f)  - R] 


[=2(t  - f 


T75 


34* 

IP" 


ds'  dt' 


34' 
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(17) 


where  R - £'  - £,  R « |£  - £'U  - (WoEq)  , 

and  n'  is  the  outward  normal  on  C (fig.  2) . The  integremds 
in  the  above  system  of  integral  equations  appear  singular  when 
R * 0 and/or  c(t  - t')  « R.  The  R ■ 0 singularity  is  only 
apparent  because  it  can  be  shown  that  R • n*  behaves  as  R^ 
where  o > 2 when  R -*■  0.  The  c(t  - t')  ■ R singularity 


-1/2 
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Is  Integrable  because  of  the  two-dimensional  integration. 

When  R X 0 and  c(t  - t')  » R simultaneously,  l.e.,  R » 0, 
t » t ' , the  factor  t - t ' provides  an  extra  zero  and  the 
R > 0 singularity  is  still  only  apparent. 

So  far  we  have  restricted  our  discussion  to  smooth 

contours.  In  this  report  we  are  Interested  in  the  numerical 

solution  for  the  problem  of  scattering  of  an  electromagnetic 

pulse  from  an  infinite  dielectric  cylinder  of  a rectangular 

cross  section.  Thus  the  behavior  of  equations  17  in  the 

vicinity  of  sharp  corners  (edges)  must  be  exeunined.  When  the 

observation  point  £ does  not  lie  at  a corner,  equations  17  are 

still  true.  This  is  so  because  it  is  well-known  that  V is  finite 

—1/2 

at  the  edge  and  B'F/Sn  varies  no  faster  than  s ' , where  s is  the 
distance  from  the  edge;  consequently  the  integrals  involving 
3'l'/3n  have  an  integrcdDle  singularity  and  are  well  behaved.  When 
the  observation  point  £ is  allowed  to  approach  a corner,  the 
factor  that  is  extracted  from  the  integral  involving  3G/3n 
is  equal  to  ± ((2/27;)  ¥ rather  than  ±1/2  ((2  is  the  interior  angle 
shown  in  fig.  3)  and  consequently  when  £ is  at  a corner 
equations  11  and  17  have  their  left-hand  sides  equal  to 
(1  - (2/27;)’?  ..  Our  numerical  treatment  for  the  pair  of 

O f 1 

integral  equations  17  will  not  allow  £ to  lie  at  a corner 
because  all  the  reference  points  are  chosen  at  the  midpoints 
of  arc  segments  as  we  will  explain  in  the  next  section.  Thus 
the  factor  to  be  extracted  is  ±1/2  and  equations  17  are  valid 
for  all  observation  points  of  interest. 

Before  we  turn  our  attention  to  the  next  section  we 
should  mention  that  in  appendices  A and  B the  validity  of 
equations  17  is  tested  analytically  by  solving  two  special 
problems  whose  solutions  are  known. 


SECTION  III 


NUMERICAL  TREATMENT  OF  INTEGRAL  EQUATIONS 

- f 

In  this  section  we  present  the  procedure  that  allows  us  to 

numerically  solve  equations  17  for  the  problem  of  scattering  of 

an  electromagnetic  wave  from  an  infinite  dielectric  cylinder  of 

a rectangular  cross  section  (fig.  4)  i.e.,  a rectangular  slab. 

The  incident  wave  is  £~polarized  in  the  x direction  and  propagates 

in  the  z direction.  (Even  though  we  focus  our  attention  on  the 

numerical  solution  for  a particular  cross  section  the  method 

we  employ  is  directly  applicable  to  other  cross  sectional 

geometries.)  In  order  to  cast  equations  17  into  a system  of 

algebraic  equations  which  we  can  solve  numerically,  we 

partition  each  of  the  four  sides  of  the  rectangle  in  figxire  4 

into  equal-sized  intervals  Ls  and  the  midpoint  of  each  line 

segment  is  chosen  as  the  reference  point  for  that  interval. 

(As  may  vary  from  side  to  side.)  In  order  to  effect  a similar 

partition  for  the  t-integration  we  observe  that  the  upper 

limit  in  equations  17  can  be  replaced  by  t since  there  can 

be  no  contribution  to  'F(£,t)  later  than  t.  Assuming  that 

the  wavefront  hits  the  front  size  of  the  rectangle  at  t = 0 

we  can  replace  the  lower  limit  of  the  t-integrations  by  zero. 

If  we  set  t = 0 in  equations  17  we  obtain  (1/2) (p,o)  = 'F^’^®(£,o) 

and  'F,  (£,o)  » 0.  There  is  no  contradiction  because  the  incident 

ij!ic 

wave  has  a smooth  wavefront  and  V (£,0)  * 0 at  all  £ on  the 
contour.  If  the  latest  time  of  interest  t is  called  T then 
we  have  a time  interval  (o,T)  that  can  be  partitioned  into 
equal-sized  Intervals  At. 

The  reference  point  for  the  time  Interval,  bounded 
by  * (j  " l)At  and  tj  « jAt,  is  tj  and  not  the  midpoint 

(j  - 1/2) At.  The  Cartesian  product  of  the  space  and  time 
partitions  produces  a lattice  of  zones.  The  reference  point 
for  the  i,j  zone  is  then  (i  - 1/2) As,  jAt,  if  for  convenience 
all  line  segments  are  of  equal  size.  Before  we  transform 
the  pair  of  Integral  equations  17  into  a system  of 
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algebraic  equations  which  can  be  solved  numerically  we 
rewrite  equations  17  in  a more  convenient  form 


r 


ti 


r t 

'i'(y»z»t)  ■ S{y,z,t) 

'o  *'C 


Zft;  y',z',t')4»(y',2',t') 


+ K,(y,2,t;  y',2',t')'l'(y',2',t')]  ds'df 


t 

'Ky^Zrt)  f (£  [K3(y»Z/t;  y' ,2' ,t')  tCy* /z' ,t' 
•b  ^ 


+ K^(y,2,t;  y',z‘,f)'l'(y',2',f)]  ds'df 
where  = 3f/3n,  ¥ = 3V3t,  S = 


(18) 


. fo  ^[Co<^  ~ ^ R] 

1 IT  r t o «-il/2 


[c^(t  - f)2  - 


''2  - -f 


®o 


t')  U[CQ(t  - f)  - R] 
[e|(t  - f)^  - 


K4  “ K2^°o  ■*"  ®i^  • 


(19) 


Our  method  of  solving  the  system  of  equations  18  is  to 
assume  that  ¥ and  * vary  so  slowly  in  space  and  time  that 
their  values  at  a point  defined  by  the  midpoints  of  the  arc 
segment  and  time  Interval  forming  a zone  provide  a good 
estimate  for  their  values  over  the  corresponding  zone.  The 
singular  nature  of  the  Icernels  forbids  us  from  malclng  the 


I 

I 

1 

I 

■j 
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seune  assun^tlon  about  them.  If  we  divide  the  circumference 
of  the  contour  C into  N segments  we  can  rewrite  equations  18 
at  t “ t j as 


- E E / ■'i  dsMf  f f Kj  d 

>'-1  “H  *k»  \t 


where  pfj  = F(p„»  nAt) , o is  the  radius  vector  to  the  midpoint 
m^  — m — m 

of  the  m^  arc  segment,  is  the  (k,l)  zone  in  the  s,t  space 
and  Kp  = Kp (£.^ » t j * f t ' ) (p  * 1,2, 3, 4). 

The  time  derivative  is  defined  as 


•i-1/2  , . 'Fj^"^)/At 

Thus  » (yJ  - 'F°)/At  - '•'J/At  since  ’1'° 

At  j » 1 equations  20  give 


/ ■'i 


ds'df  + ('pj/At)  y* 

^1 


K2  ds'df 


K3  ds'df  + (vJ/At)  f 

^1 


ds'df 


27 


i 

J 

f 


solved  to  give 

s" 


This  is  a system  of  2N  equations  with  2N  unknowns  and  can  be 

(i»l/2, . . . ,N)  in  terms  of  the  known 
quantities  St.  If  we  write  the  system  of  equations  20  at 

^ f 

j > 2,  we  again  obtain  a system  of  2N  equations  for  2N  unknowns 

T?,  in  terms  of  the  known  quantities  S?,  H'K  (i»l,2, 

* * ^ ^ i'*'  i-1/2 

Thus  we  can  march  in  time  and  solve  for  and  ' 

i'  ^i. 


for  any  i and  j in  terms  of  the  known  quantities  , ¥ 


« 


A-1/2 


(i«l,2, . . . ,N;  t ■ 1,2, . . . , j-1) . Once  we  obtain  ¥ 


j 


i-1/2 

and  ' we  can  return  to  integral  relationships  9 and  10 
and  calculate  ¥ off  the  surface  of  the  cylinder. 


Notice  that  so  far  no  restriction  has  been  placed  on 

the  relative  magnitude  between  Lt  and  R » |£'  ~ £| • Integrals 

Kp  ds'dt*  in  equation  20  represent  the  interaction  between 

the  various  spatial  segments  and  their  in^ortance  depends  on 

the  relative  magnitude  of  At  and  R (the  distance  between  points) 

as  we  will  explain  shortly.  To  make  this  clear  consider 

equations  21  written  for  j « 1.  In  this  equation  ¥?',  i.e., 

¥ evaluated  at  the  midpoint  of  the  1 line  segment  and  at 

t * At,  depends  on  4 aaid  ¥ at  the  midpoints  of  all  other 

line  segments  at  t « At.  A 2N  x 2N  matrix  has  to  be  inverted 

in  order  to  evaluate  ¥4.  However,  it  is  possible  to  choose 
1 * 4 

At  such  that  ¥^  and  in  general  ¥^  only  depends  on  ¥ and  ♦ 
evaluated  at  the  various  midpoints  at  earlier  times  without 
inverting  a 2N  x 2N  matrix,  i.e.,  each  pair  of  equations  20 
will  be  solved  for  ¥^  and  in  terms  of  S^, 

(i-1,2, . . . , j-1)  by  inverting  a 2 x 2 matrix.  The  restriction 
to  be  Imposed  on  At  is  At  £ As/2  where  As  the  size  of  line 
segments  into  which  the  sides  of  the  rectemgle  have  been 
partitioned.  (We  assume  that  all  line  segments  are  of  equal 
size;  if  not.  As  is  the  smallest  line  segment.)  The  above 
restriction  will  now  be  illustrated.  We  proceed  by  writing 
equation  21  for  i ■ 2,  idiere  the  line  segment  1 > 2 is 
depicted  in  figure  5. 
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K,  ds'dt'  + 
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K,  ds'dt' 
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K,  ds'dt' 


(22) 


Let  us  consider  in  particular  the  interaction  between  line 
segments  1 and  2.  This  interaction  is  represented  by  the 
influence  coefficient  /.  K (p,,At;  p'#t')  ds'dt'  (p»l/2,3/4). 
(Actually  K2  * » 0 if  both  the  reference  and  integration 

points  lie  on  the  same  side  of  the  rectangle  since  R * n « 0.) 
If  we  recall  equations  19  we  observe  that  contains  the 
step  function  U and  the  integration  is  to  be  performed  over 
that  portion  of  Aj^^^  that  make  U * 1,  i.e.,  c(t  - t')  - 0 

orc(t-t')  - |z-z'|^0.  This  last  inequality  defines  a 
region  of  influence  or  a light  "cone"  in  the  ct',z'  coordinate 
system  (fig.  6) . The  exact  location  of  this  light  "cone" 
depends  on  the  values  of  z and  t and  in  the  present  case 
z « As/2  and  t » At.  In  figxire  7 we  have  plotted  zone  A^^ 
(defined  by  z'  * 0/  As  amd  t « 0,  At)  for  c^At  » As/2,  As, 

3 As/2  in  a c^t',z'  space.  The  case  c 


Cj^  will  be  examined 


shortly.  Notice  that  for  c^At  > As/2  the  influence  coefficient 

dz'dt'  is  nonzero,  since  the  light  "cone"  intersects 

part  of  A.,.,  i.e.,  the  line  segment  1 influences  line  segment 

2 during  the  time  interval  At  and  consequently  V,  in  the 

1/2  ^ 

first  of  equations  22  depends  on  , i.e.,  on  4 evaluated 
at  a different  point  but  at  the  same  time  (j  « 1) . However, 
if  c^At  £ As/2  then  dz'dt' 

) since  c.At  < c.At  and  in  ec 


0 and  also  /.  K.  dz'dt' 
All  3 


0 since  c^At  < c^At  and  in  equations  22  does  not  depend 

a/2  ° * 


on  *2  (or  *2 


0) . 
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ct'  ■ ct  - (z 


■ ct  + z*  - 


Figure  6.  The  Shaded  Area  Represents  the  Light  “Cone"  or  Region  of 
Influence  for  the  Interaction  Between  Points  Lying  on  the 
Same  Side  of  the  Rectangular  Cross  Section 
(Here  the  y»0  Side) 


Next  consider  the  influence  coefficient  between  line 

segments  1 and  N,  i.e.,  line  segments  not  lying  on  the  same 

side  of  the  rectangle  (fig.  5)  and  write  equation  21  for 

i « N.  Again  w4  examine  the  influence  coefficients  /.  K 

All  P 

dz'dt'.  If  we  set  U * 1,  i.e.#  c(t  - t’)  - 0 or 

c(t  - t')  - (z'^  + ^ 0 we  find  that  the  region  of 

influence  or  light  "cone"  is  a branch  of  a hyperbola  shown 

in  figure  8.  In  figure  9 we  have  plotted  zone  (defined 
by  z - 0,  As  and  t = 0,  At)  for  c^At  * As/2,  As,  3As/2  in 
the  c^t',  z'  space.  (Notice  that  y * - As/2.)  These  plots 
exhibit  similar  features  as  those  in  figure  7,  i.e.,  the 
integral  over  A^^^^  is  zero  if  c^At  £ As/2.  In  general,  the 
presence  of  U dictates  that  the  influence  coefficent  will 
be  zero  if  t ^ t'  + R/c.  When  t * At  this  inequality  is 

satisfied  for  all  t'  (0  £ t'  £ At)  if  At  £ R/c.  The  smallest 

R is  A8/2,  since  the  reference  point  is  located  at  the  middle 
of  a line  segment,  i.e.,  cAt  £ As/2.  When  t « nAt  and  (n  - 1) 
At  £ t*  ^ nAt  we  again  obtain  the  same  criterion.  Bearing 
the  previous  discussion  in  mind,  we  can  rewrite  equations  21 
as 

’’i  - - '"I'ii 

'fl  - (^3)11  (21a) 

where  (K_).  , = /.  K ds'dt'.  Notice  that  all  influence 
P kl  A,^  p 

coefficients  are  zero  except  the  self->term8  and  (A3)j^l^< 

(The  other  two  self-terms  (K2)£]^  and  zero  because 

of  condition  R • n'«  0.)  It  is  shown  at  the  end  of  this 
section  that  the  self-terms  have  a very  simple  form,  i.e., 
(Ki)ij  > c^^t,  * c^At  Independently  of  location  and 

time.  Notice  that  the  two  Indices  correspond  to  those  in 
the  left-hand  sides  of  equations  20.  Thus  equations  21a  can 
be  solved  to  give 
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i 1 + c^/c^ 


2(^inC) 

1 


.1/2 


1 


c^At  i 


At  the  next  time  step,  i.e.,  t * 2t  equations  20  can  be 
rewritten  as 


’l  - - «l'i2  - 2 [‘"I'kl  *k''^  ^ <*2>kl 

k»l 


N 


N 


’i  - '''3>t2  ♦i''"  * E ['*3>kl  ‘J''"  " <*4'kl  ■'k''"] 

k*l 


'**“  <fl>i2  * <''3>i2  • fro"  - (fi  - 

all  fjj'^are  known  since  yO  * 0 and  'fr  are  known  from  t » t. 

^ 1/2  ^ 

Also  all  are  known  from  the  t = t step  and  the  edsove 

system  of  equations  can  be  solved  for  In  general 

by  marching  on  in  time  we  can  evaluate  in  terms 

of  (k  - 1,2,. ..,N  and  £ » l,2,...,j-l)  by 

inverting  a 2 x 2 matrix. 

The  success  of  the  above  procedure  depends  eunong  other 
factors  on  how  well  we  can  calculate  the  influence  coeffic- 
ients ds'dt'  (p  ■ 1,2, 3, 4).  Fortunately,  these 

Integrals  can  be  performed  explicitly  in  terms  of  simple 
functions.  From  equations  19  we  see  that  we  have  two  types 
of  integrals. 
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ASYMPTOTES 


Figure  10.  Geometry  Illustrating  the  Influence  of  Line  Segments  of  a Side  (here 
the  y«0  Side)  on  a Point  on  the  Same  Side.  The  Influence 
Coefficients  are  Zero  If  the  Corresponding  Zone 
Lies  Entirely  Outside  the  Light  "Cone" 
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Figure  11.  Geometry  Illustrating  the  Influence  of  Line  Segments  of  the  z«0 
Side  on  a Point  on  the  y*0  Side.  The  Influence  Coefficients 
are  Zero  If  the  Corresponding  Zone  Lies  Entirely 
Outside  the  Light  "Cone" 


f 


Figure  12.  Geometry  Illustrating  the  Influence  of  Line  Segments  of  Side  y«0 
on  a Point  on  Side  y»-d.  The  Influence  Coefficients  are 
Zero  if  the  Corresponding  Zone  Lies  Entirely 
Outside  the  Light  "Cone” 
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where  c is  or  and  A is  a zone  in  the  s,t  space.  As 

we  mentioned  ealrier,  the  presence  of  the  step  function  U 

may  alter  the  region  of  integration  since  U = 1 for  c(t  - t') 

- R _>  0 and  U * 0 for  c(t  - t')  - R < 0.  For  exeunple,  on 

the  same  side  of  the  rectangle,  say  y » 0 we  have  R » |z  - z' 1 

and  inequality  c(t-t')  - |z-z'|^0  represents  the  light 

"cone"  in  the  z',t'  coordinate  system  (fig.  6).  Thus  if  a 

zone  intersects  the  light  cone  as  shown  in  figure  10  the  area 

over  which  integrations  22  are  to  be  performed  is  the  shaded 

part  of  A.  As  we  indicated  earlier,  the  light  "cone"  does  not 

necess2u:ily  consist  of  straight  lines.  If,  for  exeunple,  the 

reference  point  lies  on  the  y >■  0 side  and  the  integration 

points  on  the  y ■ -d  or  z » 0 sides,  then  R*[(z  - z')^  + 

2 2 1/2 

or  R»[z'  + y ] respectively  and  inequality  t - t'  ^ R 

represents  a region  bounded  by  one  branch  of  a hyperbola  in 
the  z',ct'  or  y',ct'  coordinate  systems  respectively.  These 
are  the  light  "cones"  for  these  cases.  Again  if  the  zone 
intersects  the  light  "cone"  (fig.  11,  12)  only  the  shaded 
part  of  zone  A will  contribute  to  the  integration.  No  matter 
what  the  relative  position  of  the  spatial  reference  and  inte- 
gration points  is,  inequality  c(t  - t')  ^ R represents  a 
light  "cone”  or  a region  of  influence  and  in  the  appropriate 
s',t'  coordinate  system  there  are,  in  our  case,  fifteen 
possible  diagrams  for  the  relative  positions  of  the  light 


"cone"  and  a zone.  In  figiire  13  only  ten  diagreuns  are 
presented  since  intersections  of  the  light  "cone"  with  zones 
occur  symmetrically  and  only  the  right-hand  off  center  ones 
are  displayed.  'Also  the  light  "cone"  depicted  can  be  hyper- 
bolic or  straight  or  both  depending  on  the  diagram.  Notice 
that  for  each  diagram  the  integration  area  is  either  a 

rectangle  or  a triangular  region  or  a combination  of  the  two. 

Thus  we  need  only  exhibit  the  results  of  Integration  for 

diagrzuas  VII  and  XI  (fig.  14)  . First  we  recall  that  when 

the  reference  and  integration  points  lie  on  the  same  side 

A 

of  the  rect2mgle  the  inner  product  R * n*  is  zero  and  I^ 
defined  in  equation  22  is  zero.  Next,  for  convenience  but 
without  sacrificing  generality  we  choose  our  integration 
points  on  the  y • 0 side  and  the  reference  point  on  the 
z » 0 side,  i.e.,  R=  z'  e^  + y e^,  R=  (z'^  + 

(If  another  combination  is  chosen,  for  example  reference 

point  on  z s b side,  then  R«(z'-b)e  +ye, 

2 2 1/2  ~ * y 

R»((z'-b)  +y]'.  A simple  chemge  of  variables,  b * z* 

» z",  can  then  reduce  this  case  to  the  previous  one  by  appropri- 
ately changing  the  limits  of  the  z'  integration.  Similar 
agrtunents  hold  true  for  all  other  combinations.)  Referring 
to  figure  14  the  following  results  can  be  obtained. 


where  z(t')  is  obtained  by  setting  c(t  - t')  - R <■  0,  i.e., 
z(t')  - [c^(t  - t')^  - t*  - t - (1/c)  (zJ_j^  + y^)^/^ 

and 
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Similarly, 


4^ 


where  again 


f(z',t' ) 


‘k,- 


‘k-1 


'i-1 


Finally  we  calculate,  as  promised,  the  self-terms 


(Kp)ij 


c f JJ 


[c(t  - t')  - R] 

[C^(t  - t')2  - r2] 


ds 'dt ' 


(p*l,  c=c^;  p=3,  c«c^) 


where  the  reference  point  of  zone  A coincides  with  the  apex  of 
the  light  cone.  The  calculation  of  the  self-terms  does  not 
depend  on  which  side  the  line  segment  lies.  Thus  in  figure  15 
we  have  chosen  side  y » 0 and  the  above  integral  can  be 
rewritten  as 


cAt  V 


(Vu 


dudv 


o o 


(V 


2 : u2,17i 


= cAt 


where  substitutions  c(t  - t')  « v,  z 


I _ 


u have  been  made. 


In  order  to  test  the  validity  of  the  numerical  solution 
obtained  via  the  integral  equation  method,  we  calculated  'f 
at  the  middle  of  the  front  side  of  the  rectangular  slab  as  a 
function  of  time  and  also  ¥ on  the  top  side  of  the  slab  as 
a function  of  z at  a particular  Instant  and  compared  them  to 
solutions  obtained  with  the  finite  difference  method  (FDM) 
which  was  being  studied  simultaneously.  The  agreement  was 
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Figure  15.  Geometry  for  the  Calculation  of  Self-Terms 
in  the  Integral  Equation  Method 


1 


r 


excellent  and  the  relevant  graphs  (figs.  20  and  21)  are 
presented  in  section  V where  the  nximerical  results  are 
discussed.  At  that  time  we  had  decided  that  the  FDM  was 
preferedjle  for  the  calculations  of  interest  (section  VI  offers 
a comparison  of  the  two  methods)  and  we  proceeded  to  evaluate 
¥ at  points  on  the  surface  of  the  slab  using  the  FDM.  As 
additional  debugging  for  the  FDM  we  calculated  ¥ for  the 
problem  of  diffraction  by  a 90*  perfectly  conducting  wedge 
and  compared  our  results  to  the  known  exact  solution.  The 
agreement  was  again  excellent  and  this  spurred  otur  curiosity 
to  test  the  FDM  for  a wedge  with  a zero  angle,  i.e.,  a per- 
fectly conducting  half-plane.  Once  more  the  agreement  was 
excellent  and  it  suggested  the  very  interesting  possibility 
of  taclcling  open  surfaces  with  the  FDM.  In  the  next  section 
we  develop  this  method  as  applied  to  our  two-dimensional 
scattering  from  a dielectric  rectangular  slab. 


I i 


SECTION  IV 


FINITE  DIFFERENCE  METHOD 

In  order  to  apply  a finite  difference  scheme  to  the 
solution  of  our  scattering  problem  we  can  either  employ  equa- 
tions 2 and  3 for  E , H . H (option  1)  or  one  equation  for 
= 'F  (option  2),  neunely  equation  5.  Option  1 allows  us  to 
solve  for  all  field  con^nents  simultaneously  whereas  option  2 

only  produces  E and  H„,  H must  be  obtained  with  the  aid  of 
X y z 

equations  2.  However,  there  Is  a clear  advantage  of  option  2 

that  prompted  us  to  choose  it;  We  know  that  for  a perfectly 

conducting  wedge  (P.C.W. ) H and  H diverge  at  the  edge  whereas 

y * 

E^  is  zero.  We  want  to  use  the  P.C.W.  for  debugging  and  con- 
sequently we  should  use  E^  alone  which  for  our  problem  is  also 
finite  at  the  edge.  Hy  and  H^  may  or  may  not  diverge  for  a 
dielectric  body  but  in  any  case  debugging  with  the  aid  of  the 
P.C.W.  might  not  be  reliable  since  the  divergence  of  H , H 

y 2 

could  introduce  significant  errors.  We  could  try  to  provide 

special  treatment  near  the  edges  but  option  2 makes  this 

unnecessary.  In  connection  with  the  edge  behavior  one  may 

2 

wonder  whether  V E in  equation  5 diverges.  To  answer  this 

y*  * 222 

we  observe  that  this  quantity  is  equal  to  (1/c  )0  E /3t  ). 

2 ^2 

E is  finite  for  all  times  and  therefore  so  is  3 E /3t  — if 

X A X 

pathological  functions  are  excluded.  Thus  indeed 

finite.  Notice  that  as  we  cross  the  boundary,  whether  on  the 

2 

sides  or  at  the  corners,  V„_E„  suffers  a discontinuity  since 

2 *2 

c is  discontinuous  and  3 E /3t  is  continuous  (due  to  the 

* 2 2 
continuity  of  E for  all  times) . The  continuity  of  c V E 
X y ^ X 

will  allow  us,  later  on  in  this  section,  to  determine  the 
proper  finite  difference  scheme  for  this  quantity  as  we  cross 
the  boundary. 

Next  we  proceed  to  apply  the  finite  difference  method  to 
equation  5 in  a source- free  region. 


p e V 


with  appropriate  boundary  and  initial  conditions.  V stands 
for  the  two-dimensional  region  bounded  by  a contour 
(fig.  16)  emd  it  is  divided  by  the  contour  C into  an  exterior 
region  and  an  interior  region  V^.  if  V were  a homogeneous 

region  then  the  solution  of  equation  5 could  be  uniquely 
determined  at  a given  time  t and  position  £ if  the  initial 
values  of  'P  and  "P  were  known  everywhere  within  V and  'P  on 
was  known  for  all  times  up  to  t.  To  ensure  uniqueness 
of  'P(£,t)  in  our  case,  the  continuity  of  'P  and  3'P/3n  across 
the  boundary  C must  be  added  to  the  boundary  and  initial 
conditions  above.  One  may  wonder,  however,  what  this  condition 
means  when  3'P/3n  is  evaluated  at  an  edge  where  it  may  diverge. 
To  answer  this  we  recall  the  system  of  integral  equation  11 
and  observe  that  because  the  singularity  of  3'P/3n  is  integredjle 
we  can  remove  the  requirement  of  continuity  of  3'P/3n  at  the 
four  corners  (i.e.,  four  isolated  points).  We  still  obtain 
a unique  solution  for  ’P(£,t). 

Equation  5 is  a hyperbolic  equation  and  its  solution  via 
the  method  of  finite  differences  has  been  extensively  studied 
when  V is  homogeneous  (see  for  exan^le  references  2 and  3) . 

The  method  is  stedsle  and  converges  to  the  exact  solution 


Forsythe,  G.  E.  and  W.  R.  Wasow,  Finite-Difference  Methods 
for  Partial  Differential  Equations,  New  York,  John  Wiley, 
1960. 

Fox,  P.,  "The  Solution  of  Hyperbolic  Partial  Differential 
Equations  by  Difference  Methods,"  Mathematical  Methods  for 
Digital  Computers,  Edited  by  A.  Ralston  and  H.  S.  Wilf, 

New  York,  John  Wiley,  1964. 
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2 2 

when  (cT/h  ) + (cr/h  ) < 1,  where  t is  the  temporal  step 

y z 

size  and  h .h,  are  the  grid  step  sizes  in  the  y emd  z directions 

y 2 

respectively.  For  our  scattering  problem  the  finite  difference 
method  is  still  applicable  provided  the  appropriate  stability 
and  convergence  criteria  and  boundary  and  initial  conditions 
are  satisfied.  In  each  region  the  stedsllity  emd  convergence 
criteria  are:  (c^T^A^y)^  + < 1 and 

+ CCiTiAj^2^  ^ ^ outside  and  inside  the  rectangle  respectively. 

It  is  desirable  to  choose  and  h^y  = h^y,  h^^  * h^^ 

which  case  both  criteria  are  satisfied  if  ^ + 

(CoToAoz^^  < 1 since  c^  < c^.  The  boundary  and  initial 

conditions  are  those  we  mentioned  earlier  in  connection  with 

the  uniqueness  of  the  solution  of  equation  5.  (To  verify  the 

accviracy  of  the  values  at  the  interface  C we  compared  the 

fields  calculated  via  the  finite  difference  approach  and  the 

Integral  equation  method  and  the  excellent  agreement  we  obtained 

strongly  indicated  that  the  values  were  indeed  accurate.)  In 

the  present  case  the  incident  plane  wave  pulse  has  a smooth 

wavefront  and  is  assumed  to  hit  the  front  side  of  the  rectangle 

at  t » 0.  Thus  4*  and  f at  t = 0 are  known  everywhere  within  V. 

Since  derivatives  are  replaced  by  finite  differences  we  write 

y(£/0)  ■»  ['K(£fO)  - y(£f-T)]/T  and  the  initial  conditions  are 

then  translated  into  the  statement  "4'  at  t = 0,-t  is  known 

everywhere  within  V. " As  we  shall  see  later  the  finite 

difference  method  makes  explicit  use  of  this  statement.  The 

importance  of  the  boundary  condition  will  become  evident  as  | 

we  transform  equation  5 into  a system  of  difference  equations.  | 

In  order  to  accomplish  this  we  first  replace  region  V by  an  > 

orthogonal  grid  with  grid  sizes  h and  h_  as  shown  in  figure  17.  i 

y z 

Notice  that  both  emd  C coincide  with  grid  bars  and  the  points  1 

at  which  y will  be  evaluated  coincide  with  the  intersections  j 

of  the  grid  bars.  Next  we  Introduce  a temporal  step  size 
equal  to  t and  replace  * (3^y/3y^)  + 0^4'/3z^)  and 


Figure  17.  The  Grid  Used  in  the  Finite  Difference  Method. 
The  Boundaries  Coincide  with  Grid  Bars  and  the 
Observation  Points  with  the  Intersections  of 
the  Grid  Bars 
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2 2 

3 y/at  by  their  appropriate  differences  in  a region  of 
homogeneity , 

j-fCp^t)  * [’F(y  + hy,z,t)  + 4'(y  - hy,z,t) 

- 24'(y,z,t)]/h2  + 0(hy) 


(28) 

Equations  28  would  be  exact  if  4'(£,t)  were  a quadratic  function 

of  y,  z,  and  t.  In  our  subsequent  calculations  we  set  hy  » h^^ 

« h because  the  back  side  of  the  rectangle  (z  » -b)  will  be 

taken  sufficiently  far  from  the  front  side  (z  » 0)  so  it  will 

have  no  effect  on  our  field  calculations  for  the  time  periods 

of  interest.  (Thus  we  define  h by  dividing  the  front  side 

into  equal-sized  segments  and  the  top  side  is  set  equal  to 

an  integer  multiple  of  h.)  The  next  two  steps  Involve  the 

boundary  conditions  across  C and  specifying  ¥ on  The 

botuidary  conditions  across  C require  continuity  of  4'  and 

2 2 

dY/dn.  As  we  observed  earlier  in  this  section  c Vy^Y  is  also 
continuous.  These  three  conditions  will  alow  us  to  replace 
c 7^2*^  A difference  scheme  for  points  on  the  sides  of  the 
rectangle.  We  begin  with  point  £ on  side  y * 0 emd  for  con- 
venience we  choose  a coordinate  system  yz  with  origin  at  £ and 
set  t • 0.  (This  choice  for  t has  nothing  to  do  with  previous 
considerations  and  it  is  only  a matter  of  convenience.  Our 
results  will  be  valid  for  any  t.)  We  will  assume  as  we  did 


Sz' 


7(p,t) 


[7(y,z  + hj,,t)  + 7(y,z  - h2,t) 


2'P(y,z,t)]/h‘  + O(h^) 


3t‘‘ 


7(£,t)  » [7(y,z,t  + T)  + 4'(y,z,t  - t) 
- 27 (y,z,t) ]/t^  + o(t) . 


! 

4 


t 

1 

! 


I 

j. 

a 

i 

A 
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for  the  derivation  of  equations  28  that  ¥ cem  be  expanded  in 
a Maclaurin  series  about  z » 0,  y » 0. 


■ a + by  + cz  + dy^  + ez^  + fyz  + O(h^) 


y > 0 


■ a'  + b*y  + c'z  + d'y^  + e'z^  + fyz  + O(h^) 


y < 0 
(29) 


V and  3'i'/3n  « 3’l’/3y  are  continuous  as  we  cross  y * 0,  i.e., 

a*a',  b»b',  e«e',  c»c',  f » f.  To  determine  the 

2 2 

relationship  between  d and  d*  we  recall  that  c is  also 

continuous  as  we  cross  y » 0,  i.e.. 


* =i<“'  * 


(d  + e) 


In  order  to  determine  the  difference  expression  for  c 


at  y » z » ct  ■ 0 we  set 


j^V  « -ij. 

yz  ^ 


C [Af(h,0,0)  + B4'(-h,0,0)  + C'1'(0, 0,0)1 


+ lD1'(0,h,0)  + E'f(0,-h,0)  + PYIO, 0,0)1 

h^ 


- 2(d  + e) 

o 


Using  equations  29,  equation  31  gives 
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A(a  + bh  + dh  ) + B 


a - bh  + 


**  c*  * 

(d  + e)  - e h" 

L Ci  J 


+ Ca 


+ D(a  + ch  + eh^)  + E(a  - ch  + eh^) 


+ Pa  - 2(d  + e)  c^h^. 


(32) 


Setting  the  coefflcl>5nts  of  a,  b,  c,  d and  e In  equations  32 
equal  to  zero  we  obtain 

A + B + C + D + E + P-  0 

(A  - B)h  » 0 

(D  - E)h  » 0 


A + B 


2=0 


D + E + B|-j 


- 1 - 2c 


h‘  - 0 


i • e • f 


2c! 


A • B - E - D 


1 + 


8c! 


c + F ■ “ 


1 + 


(33) 
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Thus  the  difference  formula  for  c V 

o ] 

boundaries  of  the  rectangle  is 


where  y,z  is  some  point  on  the  boundary  C.  Notice  that  our 

boiindary  includes  the  four  corners  and  the  question  arises 

2 2 

as  to  what  the  difference  formula  for  c 7 7 is  at  a corner 

2 2 

The  answer  is  that  since  c 7„^7  is  continuous  we  cem  choose 

yz 

an  exterior  point  that  is  arbitrarily  close  to  the  corner 
and  use  the  first  two  of  equations  28  with  c > c . 


We  are  now  in  a position  to  replace  equation  5 by  a system 
of  difference  equations  emd  examine*  the  boundary  conditions 
for  ¥ on  the  outer  boundary  For  a point  away  from  the 

boundaries  equation  5 assumes  the  form 


where  points  (yih,z),  (y,z±h)  lie  in  the  same  medium  as 
y,z.  Notice  that  eqxiation  35  allows  the  determination  of  Y 
at  t Y-  T in  terms  of  the  value  of  Y at  £ .and  neighboring 
points  at  earlier  times.  For  a point  on  the  boundary  C 
equation  34  combined  with  the  third  of  equation  28  replaces 
equation  5.  Let  us  see  how  we  can  evaluate  Y at  a given 


point  and  time.  For  simplicity  let  us  assume  that  the  front 
side  of  the  rectangle  extends  to  infinity,  i.e. , coincides 
with  the  z > 0 plane.  This  makes  Y Independent  of  y and  the 
problem  becomes  one-dimensional  0/3y  « 0) . We  define  our 
regibn  of  interest  by  drawing  two  boundary  walls  z > z^  > 0, 
z » Z2  < 0 and  specifying  'F  on  them.  Thus  we  set  ^(Zj^^tQ)  ■ 0 
and  (Z2'^o^  * Obviously  these  conditions 

cannot  be  valid  for  any  t^  since  the  reflected  and  tremsmltted 
waves  will  eventually  reach  the  walls  z ■ Z2  and  z ■ Zj^ 
respectively.  Let  us  now  select  a point  z < 0 and  employ 
equation  35 

(cQT^/h^)t'F(z  + h,t  - t)  + 'l'(z  - h,t  - T)  - 2’i'(z,t  - t)1 

» f(z,t)  + tCz^t  - 2t)  - 2H'(z,t  - t)  (36) 

We  observe  that  depends  on  the  value  of  Y at  neighboring 

points  at  times  less  then  t.  V at  a neighboring  point  at 
t - T can  be  similarly  evaluated  by  writing  equations  analogous 
to  36  and  it  too  depends  on  the  values  of  '¥  at  neighboring 
points  at  earlier  times.  This  procedure  shows  that  VCz^t) 
depends  on  'F(z2,t  - n2T)  and  'F(Zj^,t  - nj^T)  where  |z2  - z|  ■ 
n2h  and  z^  - z » n^h.  Thus  if  the  walls  have  not  been  reached 
by  the  scattered  waves  at  t - n^T  (for  z » Z2)  and  t - nj^x 
(for  z » z^) , the  boundary  conditions  are  valid  and  so  is  the 
calculation  of  ¥ at  z and  t.  Equation  36  shows  that  ¥(z,t) 
not  only  depends  on  the  value  of  ¥ at  neighboring  points  at 
t - T but  also  on  the  value  at  ¥ at  z at  t - t and  t - 2t, 
i.e.,  a knowledge  of  these  values  is  required  which  in  turn 
depend  on  the  value  of  ¥ at  earlier  times.  This  continues 
\uitil  we  reach  the  initial  conditions.  Thus  the  difference 
scheme  works  as  follows.  First  we  write  equations  36  emd 
34  at  t • T 


{cjT^/h^)[4'(2j^,0)  + 'F(Zj^  - 2h,Q)  - 2'^{z^  - h,0)]  - 4'(2j^  - h,T) 

+ ,^(2^^  - h,-T)  - 2'f(2j^  - h,0)  2 « 2j^  - h 

I 

I 

(cjT^/h^)[H'(2  + h,0)  + 'J»(2  - h,0)  - 2'F(2,0)]  - 'i'(2,T) 

+ fCz^-T)  - 2'F(2,0)  2 > 0 

IcIt} 

-5 S-5 — 5-[y(h,0)  + 'l'(-h,0)  - 2'?(0,0)]  = 4'(0,t) 

h^d  + c‘/c‘) 

+ VCO^-t)  - 2H'(0,0)  2 » 0 

(c^T^/h^)t'P(2  + h,0)  + H'(2  - h,0)  - 24'(2,0)1  » 'F(2,t) 

+ ^(Z^-T)  - 2'F{2,0)  z < 0 

(CQT^/h^)['F(z2  + 2h,0)  + VCZj/O)  - 2^(22  + h,0)]  » + h,T) 

+ V(Z2  + h,-T)  - 2'i'(Z2  + h,0)  z » 22  + h (37) 


where  z « z.  >0  and  z * z.  < 0 are  the  two  boundary  walls 

^ * iwG 

such  that  'Hz^,t)  ■ 0,  ’F(Z2ft)  * 'f  ^2^)  • Notice  that 
equations  37  show  that  ¥(z,t)/  for  any  z (except  at  the 
boundaries),  depends  on  ¥(z  ± h,0) , ¥(2,0)  and  ¥(z,-t),  i.e., 
on  the  initial  values  (t  ■ 0,-t)  of  ¥ everywhere  within  V. 
These  values  are  ]cnown  as  we  explained  earlier,  i.e.,  equa- 
tions 37  allow  the  calculation  of  ¥(z,t)  everyvdiere.  At 
this  point  tha  boundary  conditions  are  superseded  by  the 
initial  conditions  but  they  will  manifest  themselves  in  the 
next  step  which  involves  a set  of  equations  similar  to 
equations  37  written  at  t ■ 2t  rather  than  t ■ t.  This 


59 


new  set  allows  the  calculation  of  4'(z,2t)  in  tezms  of 
¥(z  ± h,T),  ¥(Z/t)  that  were  calculated  in  the  first  step 
and  '{'(ZfO)  that  ^is  known  through  the  initial  condition. 

Notice  that  the  first  emd  last  of  equations  37  written  at 
t ■ 2t  involve  1'(Zj^,t)  and  '{'(Z2/T)  respectively,  i.e.,  the 
values  of  ¥ at  the  boundaries.  These  values  were  not  cal- 
culated in  the  first  step;  they  have  to  be  specified.  It 
is  clear  now  that  we  can  inarch  on  in  time  in  steps  of  t and 
calculate  ¥(±inh,nT)  for  any  m and  n provided  the  boundary 
conditions  are  not  violated. 

The  simple  one-dimensional  scattering  problem  above,  well 
illustrates  the  mechemics  of  the  difference  equation  method. 

Our  two-dimensional  scattering  problem  can  be  solved  similarly 
sxibject  to  appropriate  initial  and  boundary  conditions. 

The  Initial  conditions  still  require  knowledge  of  'il(y,z,0) 
and  ’F(y,z,-T)  everywhere  and  the  boundary  conditions  now 
involve  four  walls  instead  of  two,  i.e.,  y « yi»y2 
addition  to  z Z2^,Z2>  The  boundary  conditions  at  z » 
are  still  ■ 0 2md  4'(z2ft)  * (Zj^,t) . In  the  y 

direction  the  symmetry  about  the  y * yjj^  < 0 plane,  where  y^^^ 
is  the  y coordinate  of  the  middle  of  the  front  side,  allowed 
us  to  Impose  the  condition  'f  “ h,z,t)  » ^yjn  h,z,t), 
where  y * yj^  - h is  the  boundary  wall,  and  it  works  as  follows. 
The  difference  equation  at  y^,z,t  is 

Yl'?(y„  + h,z,t  - T)  + '9(y^  - h,z,t  - T)  - 2Y(y^,z,t  - x)] 

where  y depends  on  z.  Because  of  the  boundary  condition 
f(yjjj  + h,z,t  - x)  ■ T(yj|j  - h,z,t  - x) , ’•'(ynj»*»t)  can  be 
calculated  in  terms  of  the  value  of  ¥ evaluated  at  earlier 
times.  For  exao^le,  ¥(yjQ»z,2x)  is  calculated  in  terms  of 
¥(y  »z,x)  and  ¥(y«  + h,z,x)  which  were  evaluated  in  the 


first  (t  ■ t)  step  and  which  is  known  through 

t the  initial  condition.  Notice  that  the  above  boundary 

condition  is  true  for  all  times  in  contradistinction  to  the 
boundary  conditions  at  z - *j^»*2*  finally,  the  second 
boundary  condition  in  the  y direction  is  imposed  by  assuming 
. that  tCyp  + h,z,t)  ■ '•’(yp  “ h,z,t)  «diere  y ■ yp  + h is  the 

second  boundary  wall  in  the  y direction  and  y^  (as  well  as 
and  Z2)  depends  on  the  period  of  time  over  which  we  wish 
I'l  to  know  y.  This  condition  means  that  our  scattering  conflgur- 

j 

atlon  is  periodic  in  the  y direction  which  of  course  is  not 
true.  It  is  just  a convenient  condition  and  it  will  be 
violated  when  the  scattered  wave  reaches  the  wall.  Again 
the  calculation  of  f at  y#z,t  will  depend  on  the  value  of  f 
at  the  boundary  at  t - nh  where  y - y * (n  - l)h.  The 
violation  of  conditions  y(y/Z2/t)  » ‘^(y^Zj^t)  and 

y(yp+h,z,t)  ■ y (yp  “ h,z,t)  will  not  be  readily  noticed 
because  the  scattered  field  has  a smooth  wavefront  that 
builds  up  slowly  and  at  the  time  of  the  violation  the 
incident  field  may  have  a high  value.  However,  the  violation 
of  condition  y(y,Zj^,t)  « 0 corresponding  to  a perfectly  con- 
ducting wall  at  z > z^  will  be  quickly  felt  in  its  vicinity. 

We  can  easily  keep  track  of  these  violations  in  the  computer 
printout  and  discard  erroneous  results. 

As  we  explained  at  the  end  of  section  III  we  debugged  the 
finite  difference  method  by  (a)  testing  surface  field  calcula- 
tions, for  scattering  by  a dielectric  rectangular  slab, 
obtained  via  this  method  against  analogous  calculations 
obtained  via  the  integral  equation  method  and  (b)  by 
comparing  field  calculations  off  the  surface  of  a 90* 
perfectly  conducting  wedge  and  a perfectly  conducting  half- 
plane, both  illuminated  by  a plane  electrostagaetic  pulse, 
with  the  known  exact  solutions.  In  all  inatances  the  agree- 
ment was  excellent.  The  relevant  graphs  are  given  in  the 
next  section  where  all  the  ntimerical  results  are  presented. 
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Finally,  we  mention  that  the  boundary  condition  given  by 
equation  34  was  not  critical  In  applying  the  finite  difference 
scheme  to  ovir  problem.  That  Is,  by  shifting  the  grid  a little 
so  that  a grid  bar  was  just  In  front  of  the  Interface  emd  by 
applying  the  difference  scheme  given  by  equation  28  we 
obtained  results  of  coiiqpar2d)le  accuracy  to  those  obtained 
with  the  application  of  boiind£u:y  condition  34.  The  only 
difference  was  that  the  Interface  could  not  be  located  to 
within  a grid  step  size  but  this  uncertainty  cam  be  reduced 
by  making  the  step  size  smaller. 


SECTION  V 
NUMERICAL  RESULTS 

In  this  section  we  present  our  numerical  results  in  the 
I form  of  graphs  for  field  calculations  on  and  off  the  surface 
of  a rectangular  dielectric  slab  illuminated  by  plane  wave 
electromagnetic  pulses  with  smooth  wavefronts  (fig.  4) . We 
also  present  graphs,  used  for  debugging,  that  show  the 
excellent  agreement  between  calculations  obtained  via  the 
finite  difference  method  and  exact  solutions  for  diffraction 
by  0*  and  90®  perfectly  conducting  wedges.  Finally,  we 
determine  the  time  interval  over  which  our  slab  results  are 
applicable  to  the  ATLAS  I trestle  platform  problem. 

In  order  to  apply  the  finite  difference  method  to  diffrac- 
tion by  a perfectly  conducting  wedge  we  employed  the  scheme 
developed  in  the  previous  section  and  set  c^^  » 0,  i.e., 

» ».  This  is  tantaunount  to  setting  'F  * 0 for  all  times 
on  the  surface  of  the  wedge.  In  order  to  test  our  numerical 
results  we  employed  the  known  exact  solutions  for  Illumination 
by  a plane  wave  step  pulse  and  appropriately  convolved  these 
solutions  with  an  incident  pulse  of  our  choice,  i.e., 

ginc  . ^inc  , 1^  U(CQt)  U(B  - c^t) . 

® (38) 

This  pulse  starts  at  t » 0 and  terminates  at  t » B/c^  with  a 
maximum  of  100  (arbitrary  units)  at  t >■  6/2c^.  It  has  a 
smooth  wavefront  with  3'F^'^®/3t  ■ 0 at  t « 0.  A plot  of  the 
Incident  pulse  is  given  in  figure  19  with  8 >■  2.2  (arbitrary 
units) . Figures  18  and  19  exhibit  the  excellent  agreement 
between  the  numerical  results  and  the  exact  solutions.  The 
spatial  and  temporal  step  sizes  we  used  were  h « 0.05,  c^t  > 
0.025,  both  in  the  same  arbitrary  units  as  8.  Notice  that 


i 

i 
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the  observation  point  (0.1, 0.1)  is  only  two  spatial  step 
sizes  (in  both  the  y and  z directions)  away  from  the  edge 
but  the  agreement  is  still  excellent. 

The  next  two 'graphs  present  a comparison  between  the 

Integral  equation  method  and  the  finite  difference  method. 

In  figure  20  we  plotted  the  total  electric  field  at  the 

middle  of  the  front  side  of  the  slab  versus  c_t/d  where  d is 

o 

the  thickness  of  the  slab.  The  incident  field  is  still  given 

by  equation  38  with  3 - 2d  and  the  sleUo  has  a dielectric 

constemt  ® 4,  i.e.,  c^  « integral  equation 

method  we  chose  As  * d/8  and  c^At  = h/2  and  for  the  finite 

difference  method  we  chose  h = d/21  cuid  c^t  = h/2.  We  weuited 

to  make  the  finite  difference  method  as  accurate  as  possible 

in  order  to  compare  it  to  the  integral  equation  method  and 

this  is  why  we  chose  a finer  spatial  division  for  the  former 

them  the  latter  method.  Figure  21  shows  a comparison  between 

the  two  methods  for  the  total  field  evaluated  on  the  top  side 

(y  » 0)  of  the  slab  versus  z/d  at  an  instant  such  that  the 

incident  wavefront  has  just  arrived  at  z = 4d.  The  incident 

field  is  also  plotted  and  occupies  a length  of  two  sleds 

thicknesses  since  3 = 2d.  (For  this  case  we  chose  h » 2c  t » 

o 

d/21  and  As  = 2c^At  = d/6.)  The  agreement  between  the  integral 
equation  method  and  the  finite  difference  scheme  is  excellent 
and  this  served  as  debugging  for  both  approaches. 

The  subsequent  graphs  present  field  calculations  via  the 
finite  difference  method  at  points  off  the  surface  of  the 
slab.  It  was  determined  that  for  these  calculations  the  finite 
difference  method  was  superior  to  the  integral  equation  method 
because  of  a lesser  use  of  computational  resources  (see  sec^ 
tlon  VI) . The  Incident  field  was  chosen  as  a fast  rising 
and  slowly  decaying  pulse  in  order  to  emulate  the  electromag- 
netic pulse  arising  from  an  exoatmospheric  nuclear  detonation. 
Such  a pulse  is  appropriate  for  studying  the  effect  of  the 


TOTAL  FIELD  AT  THE  MIDDLE 
OF  THE  FRONT  SIDE  OF  SLAB 


Strength  are  Arbitrary 


i 

I 

I 

I 


i 


1 


ATLAS  I simulator  trestle  platform  on  the  Incident  electromagnetic 
pulse.  We  will  discuss  the  applicability  of  our  results  to  the 
ATLAS  I trestle  platform  problem  later  on  in  this  section.  The 
incident  pulse  is 'given  by 

ginc  ^ ,^inc  _ 331361(d/c^t)^  ^-121. 5 (d/cot)  . (39) 

This  pulse  has  a smooth  wavefront  (all  derivatives  vanish  at  t ■ 0) 
and  it  reaches  a maximum  of  101.422  (arbitrary  units)  at  c^t  " 9d. 
However,  at  c^t  ■ 3d  the  pulse  has  only  risen  to  0.017.  For  this 
reason  we  arbitrarily  set  t ■ 0 when  the  field  strength  at  z ■ 0 
is  0.017.  Thus  the  effective  rise  time  is  Sd/c^,  i.e.,  from  0.017 
to  peak.  If  we  adopt  the  definitions  given  in  references  4 and  5 
then  for  our  pulse  t^„  % 3.03  d/c^  end 

^10-90  " rise  time 2.98  d/c^.  The  reason  we  chose 

equation  39  for  our  incident  pulse  is  that  the  usual  double  expon> 
entlal  employed  by  other  workers  to  represent  the  true  EMP  has  a 
discontinuous  first  derivative  at  t « 0.  (Notice  that  we  could 
also  have  employed  the  Inverse  double  exponential  waveform 
^inc  ^ ’l'^[ (exp(-a(t-tQ))  + exp(S(t-t^))  ]”^  (refs.  4 and  5)  for  which 

all  derivatives  at  t - 0 exist.  A suitable  choice  for  t^  should 

inc  ° 

make  ¥ at  t * 0 as  small  as  desired.) 

We  have  plotted  the  total  electric  field  at  points  above  the 
slab  versus  the  distance  from  the  origin,  with  the  vertical 
distance  from  the  top  of  the  slab  and  the  dielectric  permittivity 
as  parameters.  It  is  assximed  that  the  effective  wavefront  (i.e., 
a field  strength  of  0.017)  hits  the  edge  at  t - 0.  In  figure  22 
the  graph  for  z ■ 5d  starts  at  c^t/d  6 and  at  this  instant  the 
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5.  Castillo.  J.  P.,  K.  C.  Chenand  and  C.  E.  Baxan,  Relation  of  Rise  j i 

Time  Definitions  for  Various  Waveforms.  Atlas  Memo  20,  Air 
Force  Weapons  Laboratory,  November  1976. 
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Figure  23.  Incident  and  Total  Electric  Field  (In  Arbitrary  Units)  Versus 
Normalized  Time  at  z ■ 5d  with  e.|  '>8eg.  The  Normalized  Distance 
Above  the  Slab  Is  the  Parameter 
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Figure  24.  Incident  and  Total  Electric  Field  (In  Arbitrary  Units)  Versus 
Normalized  Time  at  z ■ 15d  with  ‘^e..  The  Normalized 
Distance  Above  the  Slab  Is  the  ^Parameter 
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Figure  25.  Incident  and  Total  Electric  Field  (In  Arbitrary  Units)  Versus 
Normalized  Time  at  z ■ 15d  with  ci«  8e..  The  Normalized 
Distance  Above  the  Slab  Is  the  ^Parameter 
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Figure  27.  Incident  and  Total  Electric  Field  (In  Arbitrary  Units)  Versus 
Normalized  Time  at  z • 25d  with  ei*8e..  The  Nonnallzed 
Distance  Above  the  Slab  Is  tne  Parameter 


incident  field  has  a value  of  2.6,  i.e.,  the  effective  wavefront 
has  already  passed  the  point  z ■■  5d  for  a length  of  time  At  ■ d/c^. 
Graphs  at  points  further  away  from  the  edge  also  show  the  actual 
times  corresponding  to  the  incident  wave  sweeping  by.  Thus  at 
z - 15d  we  start  at  t ■ 16d/Cg  and  at  z - 25d  at  t ■ 26d/c^.  To 
obtain  the  results  displayed  on  the  graph  through  z - 25d  we 
pushed  the  memory  requirements  under  GDC  7600  FTN4  to  the  limit 
and  therefore  we  were  forced  to  use  a relatively  coarse  grid  size, 
h * d/7.  Thus  we  cannot  guarantee  the  accuracy  of  the  results. 
However,  results  for  h - d/21  (at  z » 5d)  and  h ■ d/ll  (at  z ■ 15d) 
were  similar  to  those  at  h * d/7  over  the  common  region  in  space 
and  time,  the  maximum  deviation  being  1 to  2 percent  of  the  peak 
incident  field. 

The  plots  show  that  the  presence  of  the  slab  can  distort  the 
incident  field  significantly.  On  the  top  surface  of  the  slab  or 
very  close  to  it,  the  total  field  reaches  a maximum  that  is 
shifted  in  time  relative  to  the  maximum  of  the  incident  field. 

This  maximtim  is  also  larger  than  the  peak  value  of  the  Incident 
field.  The  behavior  of  the  total  field  can  be  qualitatively 
wderstood  if  we  take  into  account  the  secondary  wave  within  the 
slab  propagating  with  a speed  c * c^.  The  larger  the  dielectric 
permittivity  the  slower  the  secondary  wave  is  (i.e.,  smaller  c^) 
and  the  total  field  reaches  its  peak  later.  This  can  be  seen  from 
the  graphs  for  and  • Se^.  As  the  observation  point 

moves  upward  the  total  field  tends  to  exhibit  two  humps  until  it 
is  sufficiently  high  where,  due  to  the  diminishing  influence  of 
the  slab,  resembles  the  incident  field.  As  the  observation  point 
rescinds  from  the  edge  the  presence  of  the  slab  becomes  more 
pronounced  and  one  must  reach  progressively  higher  observation 
points  (larger  y/d)  before  the  influence  of  the  slab  has  diminished. 
Thus  at  z * 5d,  y ■ 6.5d  (e^  * field  exhibits,  in  some 

approximate  manner,  the  same  distortion  as  the  field  at  z * 25d, 
y - 16. 5d  (Ej^  - 8e^)  . 


1.  THE  ATLAS  I TRESTLE  PLATFORM 


Our  Cwo- dimensional  scattering  from  a rectangular  dielectric 

slab  can  serve  as  a model  for  studying  the  Influence  of  the  wood 

support  structure  (trestle),  of  the  ATLAS  I simulator,  on  the 

waveform  of  the  simulated  EMP.  We  will  only  focus  our  attention 

on  the  platform,  l.e..  Ignore  the  rest  of  the  support  structure. 

(See  references  6 and  7 for  a study  of  the  reflection  of  a plane 

wave  from  the  rest  of  the  support  structure.)  This  platform  Is 

depicted  In  figure  28.  First  notice  that  the  effective  rise  time 

of  the  Incident  pulse  given  by  equation  39  Is  6d  and  since  the 

platform  thickness  Is  approximately  2 feet  this  effective  rise 

time  is  translated  into  12  nsec.  Whereas  t „ and  t.^  q-  are 

max  rise  10-90 

approximately  equal  to  6 nsec.  Thus  our  pulse  Is  faster  rising 
(also  faster  decaying)  than  the  actual  pulse  to  be  fired  In  the 
ATLAS  I simulator  (see  ref.  8 for  the  waveforms  obtained  In  the 
pulser  test  fixture  (PTF)  with  the  ATLAS  prototype  pulser  module 
and  ref.  9 for  a summary  of  the  final  results  of  all  testing 
performed  on  the  ATLAS  prototype  pulser  module  In  PTF  and  also  the 
Influence  of  additional  diffraction  and  relfectlon  effects  on  the 
waveform  of  the  pulse . ) This  difference  in  the  rise  time  and  other 
pulse  characteristics  between  our  pulse  and  the  one  to  be  fired 
Into  the  working  volume  of  ATLAS  I makes  our  quantitative  results 
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not  directly  applicable;  our  results,  however,  can  still  provide 
information  on  the  influence  of  the  platform  on  the  waveform  of  a 
fast  rising  and  slowly  decaying  pulse.  Next  we  observe  that  our 

f 

two-dimensional  study  is  valid  as  long  as  it  is  applied  to  observa- 
tion points  that  have  not  been  reached  by  the  diffracted  waves  due 
to  the  edges  P',  P".  Choosing  our  observation  points  on  the  x • 0 
plane  the  maximvim  or  clear  time  Interval  over  which  our  two- 
dimensional  model  is  valid  is  the  time  that  elapses  from  the  instant 
the  incident  field  reaches  the  observation  point  P until  the  instant 
the  diffracted  fields  from  edges  P' , P"  reach  P (fig.  28),  i.e.. 
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Thus,  assuming  that  for  the  wood  platform,  c^t/d  in 

figure  22  extends  throughout  the  indicated  region,  in  figure  24 
approximately  throughout  the  indicated  region  and  in  figure  26  up 
to  an  average  of  26+17-43.  In  any  case  the  graphs  show  che 
distortion  of  the  incident  waveform  due  to  the  presence  of  the 
platform,  over  a large  portion  of  the  incident  pulse,  well  past 
its  rise  time. 


9.  Baum,  C.  E.,  D.  E.  Higgins  and  D.  V.  Giri,  Pulser  Test  Results 
and  Preliminary  Estimation  of  Transient  Electric  Field  Waveforms 
in  ATLAS  I.  ATLAS  Memo  18,  Air  Force  Weapons  Laboratory,  Oct.  1976. 
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SECTION  VI 


COMPARISON  OF  THE  INTEGRAL  EQUATION  AND  FINITE 
DIFFERENCE  METHODS 


In  this  section  a comparison  is  drawn  between  the  integral 
equation  method  and  the  finite  difference  scheme  based  on  the 
required  computational  resources  and  degree  of  overall  sim- 
plicity. 

1.  MEMORY  REQUIREMENTS 

a.  Finite  Difference  Scheme 

In  order  to  ensure  that  the  imposed  boundary  conditions 

discussed  in  section  IV  do  not  affect  the  accuracy  of  the 

field  strengths  con^)uted  in  the  region  outside  the  dielectric 

body,  it  is  necessary  that  the  grid  extend  a distance  c^T  in 

the  positive  and  negative  z-directions  and  in  the  distance  above 

the  platform,  where  T is  the  length  of  the  time  from  when  the 

incident  field  first  hit  the  body  until  the  latest  time  of 

interest.  (Due  to  symmetry  about  the  plane  y * < 0, 

where  y„  is  the  y coordinate  of  the  middle  of  the  front 
m 

side,  only  the  y > y^^^  region  need  be  considered.)  If  the 
only  fields  of  interest  are  those  on  or  inside  the  body,  the 
grid  need  only  be  extended  to  a distance  c^T/2  instead  of  c^T. 
(This  can  be  understood  by  recalling  the  influence  of  the 

boundary  conditions  discussed  in  section  IV. ) Thus  the  grid 

. 2 
size  must  be  2c^T(c^T  + d/2)/h  in  the  former  case  and 

CqT(c^T/2  +d/2)/h^  in  the  latter  case  where  h is  the  grid 

step  size.  The  finite  difference  algorithm  requires  knowledge 

of  the  fields  at  the  present  time  step  and  the  previous  time 

step  in  order  to  confute  the  fields  at  the  next  time  step. 

However,  the  field  strength  at  a given  grid  point  at  the 

previous  time  step  is  only  needed  to  compute  the  value  at 

the  next  time  step  at  the  same  grid  point  and  consequently 

1 t is  only  necessary  to  provide  two  storage  locations  for 

• srid  point. 


b.  Integral  Equation  Method 

The  Integral  equation  method  permits  a tradeoff  in  costs. 

Instead  of  recomputing  the  influence  coefficients,  discussed 

in  section  III,  for  each  time  step,  it  is  possible  to  store 

them  for  future  use  and  compute  only  those  coefficients 

relating  the  present  time  step  to  the  first  time  step.  We 

can  do  this  because  the  kernel  of  the  Integrals  defining  the 

influence  coefficients  depends  on  t - t*  and  not  t or  t' 

individually.  Thus  there  is  a tradeoff  in  that  presumcdsly 

memory  references  are  faster  than  evaluation  of  transcendental 

functions  (resulting  from  the  explicit  calculation  of  the 

double  Integrals  defining  the  influence  coefficients)  and 

the  logic  required  to  determine  whether  or  not  a zone  can 

influence  the  reference  point.  Under  the  assumption  that 

this  tradeoff  of  memory  for  computer  time  will  be  used, 

2 

4N  (T/At)  memory  location  are  required  to  store  the  influence 

coefficients  where  N is  the  number  of  line  segments  and  At 

is  the  time  step  size.  This  number  of  memory  locations  can 

be  reduced  by  Implementing  existing  symmetry  relations  (e.g., 

interaction  of  neighboring  line  segments  is  independent  of 

location  provided  that  the  segments  lie  on  the  same  side) 

but  the  basic  dependence  remains  unchanged.  It  should  be 

noted  here  that,  for  the  angle  of  incidence  of  the  incident 

field  we  are  considering  in  this  report,  if  the  length  of 

the  top  and  bottom  sides  of  the  rectangle  is  greater  than 

c^T,  N can  be  reduced  to  Ignore  the  region  of  the  rectangle 

that  will  not  be  illuminated  by  the  incident  field  by  time 

c^T.  Thus  for  such  cases  the  number  of  memory  locations  for 

the  Influence  coefficients  depends  on  three  additive  terms 

2 3 

proportional  to  T,  T and  T respectively,  and  if  the  front 

3 

side  is  much  smaller  than  c_T,  T dominates.  In  addition 

o , 

the  integral  equation  method  requires  knowledge  of  y and  4 
for  all  previous  times,  l.e.,  2N(T/At)  memory  locations 
whether  or  not  the  influence  coefficients  are  stored. 
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2.  COMPUTATION  TIME 

a.  Finite  Difference  Scheme 

Ignoring  th^  boundary  points  which  require  special 

treatment  each  grid  point  requires  the  same  length  of  time 

at  each  step.  From  the  form  of  the  difference  scheme  which 

replaces  the  differential  equation  at  each  grid  point  we 

can  see  that  three  multiplications,  six  additions  and  eight 

array  references  are  needed  for  the  calculation  of  the  field 

strength  at  given  grid  point  at  the  next  time  step.  These 

operations  must  be  performed  for  each  grid  point  for  each 

2 2 

time  step,  i.e.,  c^T  (c^T  + d/2)Tj/h  overall  computation 
time  is  required  over  the  time  of  interest  T,  where  TJ  is 
the  time  length  of  the  operations  indicated  above.  If 
c^T  >>  d/2  the  time  required  is  proportional  to  T . This 
result  may  be  reduced  by  factor  of  three  by  recognizing 
that  in  calculating  the  field  strength  at  a grid  point  at 
t the  finite  difference  technique  need  not  be  applied  to 
points  at  a distance  greater  than  c^t  from  the  body  because 
the  scattered  field  has  only  traveled  a distance  c^t  and 
the  field  strength  at  grid  points  outside  this  region  is 
that  of  the  incident  field.  If  the  total  period  of  interest 
is  T then  the  computation  time  is  approximately  c‘T  T*/3h 

O M 

if  c_T  >>  d/2.  (The  factor  of  three  comes  from  t dt  » 

o o o 

TV3.) 

b.  Integral  Equation  Method 

As  we  explained  in  our  discussion  of  the  computation 
time  for  the  finite  difference  scheme  each  grid  point  requires 
the  seune  length  of  time  at  each  time  step  independently  of 
the  parameters  of  the  system,  i.e.,  length  of  time  the  body 
has  been  exposed  to  the  incident  field,  size  of  body  and 
step  size.  Unlike  the  finite  difference  scheme  the  central 
processor  time  for  the  Integral  equation  method  depends  on 
all  those  parameters.  Asstuning  now  that  the  influence 
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coefficients  have  already  been  calculated  and  stored  the 

calculation  of  the  fields  at  time  t requires  approximately 
2 

4N  (t/At)  multiplcations  and  additions.  Thus  to  extend  the 

2 2 

calculation  to  time  T requires  2N  (T/At)  multiplications 

and  additions,  and  if  Ti  is  the  time  of  one  multiplication 

^ 2 2 
plus  one  addition  the  computation  time  is  2N  (T/At)  T^. 

The  calculation  of  the  influence  coefficients  requires  the 

2 

calculation  of  4N  (T/At)  double  integrals.  Even  though  for 

our  geometry  these  integrals  can  be  done  analytically,  the 

necessary  evaluation  of  the  resulting  transcendental  functions 

2 2 

will  contribute  a significant  amount  to  2N  (T/At)  T^  unless 
T/At  is  large.  It  should  be  noted  that  the  edsove  discussion 
assumed  that  only  the  surface  fields  were  of  interest.  For 
points  off  the  surface  of  the  body  the  calculation  of  ¥ au 
a given  time  t requires  2N(t  - R/c)/At  additional  memory 
locations  for  the  influence  coefficients  and  2N(t  - R/c)/At 
additional  multiplications  and  additions  where  R the  average 
distance  of  the  reference  points  on  the  surface  of  the  body 
to  the  observation  point  off  the  surface.  (If  'F  at  the 
same  point  is  desired  for  a series  of  equally  spaced  time 
steps,  the  required  number  of  Integral  evaluations  will 
depend  only  on  the  greatest  time  of  Interest,  not  on 
the  number  of  values  wemted.) 

3.  CONCLUSIONS 

As  we  mentioned  in  the  beginning  of  this  section  the 
comparison  between  the  two  methods  will  be  based  on  the 
degree  of  overall  sinpllcity  etnd  the  required  computational 
resources.  With  respect  to  deriving  the  necessary  equations 
and  casting  them  into  a form  suitable  for  numerical  calcula- 
tion the  finite  difference  scheme  is  much  sin^ler  than  the 
Integral  equation  method  and  this  simplicity  is  definitely 
a great  advantage.  Both  methods  can  be  made  comparable 
in  accuracy  but  the  required  computational  resources  for 
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the  two  methods  depend  differently  on  the  length  of  time 

of  interest  and  this  factor  influences  the  regions  over 

which  one  method  can  be  superior  to  the  other.  To 

understand  this'^int  recall  that  the  finite  difference 

2 2 2 

method  requires  approximately  4c  T /h  memory  locations 
. 2 3*2  ^ 

1 and  c^T  Tj^/3h  compuatation  time  whereas  the  integral 

equation  method  requires  approximately  4N^ (T/At)  (or 
3 3 

T /(At)  if  the  incident  wave  has  illuminated  a length 
c^T  of  the  body  and  c^T  >>  d/2)  memory  locations  and 

2N^(T/At)^T^  computation  time.  Thus  assuming  that  2c  At  h 

^ 2 ° 

we  understand  that  when  T/At  < N the  finite  difference 

scheme  is  preferable.  However,  for  large  times,  i.e., 

2 

T/At  > N the  integral  equation  method  requires  fewer 
computational  resources  and  should  be  considered  if  the 
simplicity  of  the  difference  equation  scheme  is  overridden 
by  an  appreciably  smaller  cost  for  the  integral  equation 
method.  If  the  field  strength  must  be  calculated  at  points 
far  away  from  the  body,  the  integral  equation  method  is  more 
practical  since  the  finite  difference  scheme  requires  propa- 
gating the  field  from  the  body  to  the  region  of  interest 
whereas  for  the  integral  equation  method  the  distance  from 
the  body  has  no  effect  on  the  required  resources. 

Our  conclusions  about  computer  resources  were  based  on 
the  two-dimensional  model  of  the  wave  equation.  For  a 
three-dimensional  model  the  finite  difference  method  would 
again  be  relatively  insensitive  to  the  dimensions  of  the 
body.  Memory  requirements  would  vary  as  T hence  total 

4 

time  will  vary  as  T . For  the  integral  equation  method 
it  is  no  longer  necessary  to  store  Influence  coefficients 
and  field  strengths  for  all  time  but  only  for  time 

steps  where  r^  is  the  maximum  diameter  of  the  body.  The 
number  of  patches  will  vary  as  L^/h^  where  L and  h are 
characteristic  dimensions  of  the  body  and  spatial  step  size. 
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hence  the  matrix  of  Influence  coefficients  requires  of 
2 2 2 

the  order  of  (L  /h  ) (r  /c^At)  memory  locations.  Similarly 

^ 2 2 

the  number  of  field  strengths  required  is  (L  /h  )(r  /c^At). 

2 2 nr  1 

The  time  per  partch  per  step  varies  as  (L  /h  ) (r^^/c^Lt)  hence 
total  time  varies  as  (L^/h^)  (rjjj/Cj^At)T.  From  the  above 
brief  discussion  we  understand  that  we  can  draw  similar 
conclusions  about  the  relative  merits  of  the  two  methods 
as  for  the  two-dimensional  con^arlson. 
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APPENDIX  A 

REFLECTION  FROM  A DIELECTRIC  HALF  SPACE 

r 

We  assume  that  the  incident  electric  field  has  the  form 
ginc  ^ ^inc  4'^“°  » f(t  - y^/c^)  - ftt  - l/c^Cy  sin  a - 

z cos  a)]  (fig.  Al) . Thus  at  the  boundary  surface  (z  * 0) 
inc 

» f(t  - y/c^  sin  a).  It  is  reasonable  to  assume  that 
the  total  surface  field  ¥ and  dV/dn'  will  only  depend  on 
t - y/c^  sin  a,  i.e.,  as  the  incident  surface  field  propagates 
to  the  right  the  total  surface  field  keeps  pace  with  it.  By 
noting  that  (£'  - £)  • n **  (y'  - y)ey  • » 0,  equations  17  can 

be  written  as  follows 


(A-1) 


Invoking  causality  we  understand  that  the  upper  limit  of  the 
t*  integration  can  be  replaced  by  t.  (The  lower  limit  is  of 
no  interest  for  the  calculation  in  this  appendix.)  Setting 
D ■ 1 we  obtain  y'  ■ ±c(t  - t')  as  the  limits  of  the  y'  inte- 
gration . 

Next  we  set  t - t*  ■ t and  the  limits  of  the  t and  y' 
integrations  are  (o,*)  and  (-ct,ct)  respectively.  Recalling 
that  aVSn  only  depends  on  t'  - (y'/Cg)  sin  o we  make  the 
following  orthogonal  transformation 


u 


H 


V' 

T + sin  o 
°o 

-T  sin  0+5^ 


(A-2) 


Thus  34'/3n  is  a function  of  t'  -(y'/Cq)  sin  a « t - t - (y'/c^) 
sin  a,  i.e./  a fvinction  of  t - u which  is  independent  of  v 
since  t is  just  a constemt.  This  will  allow  us  to  do  the 
V- integration  explicitly.  Noting  that  dudv  « l/c^d  + sin^  a) 
drdy*  and  referring  to  fig\ire  A2  we  can  rewrite  equations  A-1 


'J'(o,t) 


'F^*'°(o,t) 


1 f dv 


’f  (o,t) 


dv ^ 

(AjV^  + BjV  + C,)' 


(A- 3) 


where  L. 


It  sin  a 
sin  a ± 1 

Y T sin  g 
Y sin  a ± ] 


■ -cos  a 


■ -4  sin  o 

■ cos^  a 


Y^  sin^  0-1 

-2(1  + Y^)  sin  a 

2 2 
Y “ .sin  o. 


If  we  perform  the  v-integration  we  obtain 


|''i'(o,t)  - 4'^“®(o,t)  - 


cos  a 


I f(o,t)  - J 


(l  - sin^ 


(A-4) 


where  I 


Solving  Equation  A-4  for  we  obtain 


'f(o,t)  ^ 2y  cos  g 

77oTa  + (l  - ,i„2  a)  ■ 


Y » ^ (A-5) 


i.e.,  the  well-known  result  for  the  transmission  coefficient. 
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Mh* 


Mh* 


'F(t,z=-d) 


A(t)  = 2f(t)  - c^B(t)  - c^D(t  - d/c^)  - C(t  - d/c^) 

A(t)  « c^B(t)  + c^D(t  - d/c^)  + C(t  - d/c^) 

C(t)  » 2f(t  - d/c^)  - c^B(t  - d/c^)  - c^D(t)  - A(t  - d/c^) 

C(t)  « c^B(t  - d/c^)  + C^D(t)  + A(t  - d/c^)  (B-2) 
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In  order  to  solve  Equations  B-2  we  trzuisform  them  into  the 
frequency  domain: 

ikod  ilCo^ 

A((d)  + B(w)  +'C(aj)  e - D(u)  e ■ 2f(w) 


1 


ikid 


Aim)  - c^BCw)  - C(a»)  e + c^D((t))  e 


ik^d 


ik  d ik^d 

A((t>)  e ° + c B(«)  e ° + C(u)  - c D((o)  ■ 2f((i))  e 


ik.d  ikjd 

-A{a))  e ^ - c^B(u)  e ^ + C(u)  + c^D((i))  - 0 (B-3) 

Solving  Equations  B-3  for  ACu)  we  obtain  the  well-known 
result  for  the  reflection  coefficient  R((d) 


R((i))  ■ A(u)  - f («) 


(1  - Y^)  i sin  g 


2y  cos  o - i(l  + Y ) si*'  oi 


idiere  y * ®i/°o  “ ^o-^i* 
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